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I. INTRODUCTION 


System identification technology has been used successfully 
for many vehicles. Because of their large number of degrees of 
freedom and complex aerodynamic interactions, rotorcraft have 
always presented a special challenge to system identification 
methods. A completely integrated methodology has been developed 
under this NASA contract to solve this difficult problem. This 
methodology has also been translated into a user-oriented series 
of computer programs. This volume provides basic guidelines for 
efficient and effective use of one of these computer programs. 

Figure 1.1 shows a schematic flowchart of the overall data 
processing technique for rotorcraft. The first step in this pro- 
cedure is state estimation and instrument calibration. This is 
implemented by the computer program DEKFIS (Discrete Extended 
Kalman Filter and Smoother) which implements an extended Kalman 
filter/smoother using the Friedland-Duffy formulation. Instrument 
biases and scale factors are estimated at this stage, together with 
any state which is not measured directly. The second step involves 
estimation of the mathematical model of various forces, moments and 
interchanges. This is implemented in OSR (Optimal Subset Regres- 
sion) computer program which uses a regression technique. Accurate 
estimates of parameters are obtained in the final step. One of 
two computer programs is used for this purpose. SCIDNT implements 
the maximum likelihood method for linear systems and NLSCIDNT ex- 
tends the method to nonlinear rotorcraft models. 

Accuracy of parameter estimates may be improved by using 
flight test inputs based on the input design program, INDES. 

DEKFIS is a discrete extended Kalman filter/smoother program 
formulated for aircraft and helicopter state estimation and data 
consistency. DEKFIS is currently set up tv pre-process raw test 
data by removing biases, correcting scale factor errors and 

The contract research effort which led to the results in this report was 
financially supported by the Structures Laboratory, USARTL (AVRADCOM), 

NASA Langley Research Center and NASA Ames Research Center. 




PROVtOING INPUT TIME HISTORIES 
WHICH OPTIMIZE INFORMATION ON 
DESIRED VEHICLE PARAMETERS 


VEHICLE/ROTOR 
STATE ANO OUST 
ESTIMATES 
III II ■ 


INITIAL SELECTION 
OF MODEL OOF: 
FUSEUCE/FLAP/LAG 
VARIABLES. ETC. 


PROVIDES FILTERED STATES 
PROVIDES QUANTITATIVE DATA 
CONSISTENCY EVALUATION 
PROVIDES OUST TIME HISTORY 
(IF GUSTS PRESENT) 


MODEL S' 
OETERM 

1 

TPUC7URE 
I NATION 

o: 

1 

SR 

r~* ~ ■* 


SELECTS REQUIRED DO* 
SELECTS REQUIRED STABILITY 
ANO CONTROL DERIVATIVES 
AND RANKS THEIR IMPORTNACE 
PROVIDES START-UP VALUES 
FOR ML PROGRAM 


LlhtAR MAXIMUM 
LIKaiHOOO 
PARAMETER 
IDENTIFICATION 


SCIONT 


NONLINEAR MAXI- 
MUM IIKELIHOOO 
PARAMETER 
IDENTIFICATION 


NLSCIDNT 


f: 


PROVIOES UNBIASED. 
EFFICIENT PARAMETER 
ESTIMATES 

9ROVIOES 3 3AROS ON 
PARAMETER ESTIMATES 
PROVIOES FURTHER 
INFORMATION ON MOOEL 
STRUCTURE USING 
INFORPWTION MATRIX 


Figure 1.1 


Integrated Rotorcraft System 
Identification Procedure 


2 




providing consistency wit* the aircraft inertial and kinematic 
equations. Hence, in its present formulation, DEKFIS will tend to 
suppress errors introduced by the measurement system and its dy- 
namics, without affecting the overall system dynamics. Because of 
this, it is possible to separate the program into three separate 
estimators: the fuselage/gust estimator, the rotor state estimator 

and the RSRA estimator. This not only reduces the overall dimen- 
sionality and computational burden, but provides additional flexi- 
bility as shown in Table 1.1. 

Table 1.1 
Options 



FUSELAGE/GUST 

ESTIMATOR 

ROTOR STATE 
ESTIMATOR 

RSRA 

ESTIMATOR 

FIXED MING 
AIRCRAFT 

/ 



SINGLE ROTOR 
HELCIOPTER 

/ 

/ 


MULTI ROTOR 
HELICOPTER 

/ 

Repeat for 
each rotor 


RSRA- FIXED 
WING CONFIG. 

/ 


/ 

¥ 

RSRA-COMPOUND 

CONFIG. 

V 

/ 

/ 

WHIRL STAND 
DATA 


✓ 



The overall estimator has the same basic estimation algorithm 
for all the estimator blocks- -only the subroutines containing the 
state and measurement equations and their partial derivatives 
change. The algorithm is a Friedland-Duffy extended Kalman filter 
with a locally iterated smoother for the filter pass, and a 
fixed interval smoother to provide gusts, process noise and 
initial conditions when desired. 
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II. PROGRAM DESCRIPTION 


2 . 1 PROGRAM STRUCTURE 


The overall logic of the DEKFIS program is presented in Fig- 
ure 2.1. It is intended to give the user a general overview of 
the program. The subroutine calling structure is given in Figure 
2.2. The user will find this figure useful in further under- 
standing the program flow and in constructing an overlay structure 
should one be needed. 

The functions of the more important subroutines are described 
below. 

DEKFIS is the main routine. It performs no calculations, 
but calls four subroutines- - INPUT, DIRCTR, FIS, and 
OUTPUT- -which accomplish the computational and input/ 
output tasks. DEKFIS sets up the maximum dimensions 
of the major arrays, and passes the appropriate STATE 
and MEAS subroutine addresses to DIRCTR and FIS, depend- 
ing on whether the fuselage/gust estimator, the rotor 
state estimator, or the RSRA state estimator is being 
executed. 

BLOCKD is a block data routine. Its major function is to 
initialize arrays which contain labels for the primary 
states, secondary states, measurements and noise sources 
for input/output purposes. 

INPUT handles all the program's card and mass storage input, 
checks for errors in the input, initializes various 
flag arrays, and prints the title page and a summary 
of this information. 

INREAD is called by input to read the time histories of the 
measurement and control variables. 

PRTDAT prints out the measurement and control time histories. 

SETDIM sets the dimensions and pointers for the component 

matrices of the D and R arrays. The dimensions are for 
the compressed matrices and vectors for the problem 
currently being executed. 
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2 Subroutine Calling Structure 



DIRCTR performs the main computational tasks for the discrete 
extended Kalman filter and locally iterated smoother 
portion of the program. Its major features are a bias 
free filter, a bias free local smoother, a bias filter, 
a bias correction subroutine and a composite state local 
smoother. DIRCTR performs the initialization of the 
STATE, MEAS and other routines. It then sequences 
through all the data points yielding the filter estimate 
at each point. 

BFFLTR performs the bias free filter calculations. It is 
called at each time point. The bias free filter is 
simply an extended Kalman filter for the primary states 
alone (i.e., the secondary or bias states are held 
constant). Hence, it performs prediction and update 
calculations for the primary states and the associated 
covariance. 

IN0VAT is called by BFFLTR to calculate the innovations for 
a given time point. 

STATRN is called by BFFLTR to calculate the state transition 
matrix for a given time point. 

INTEGX is called by BFFLTR to integrate the state equation 

to the next time point for the primary state prediction. 

STATE contains the itate equations for the appropriate esti- 
mator. It has entry points XDOTV, FMAT, GMAT, GAMAT 
and BMAT. XDOTV contains the nonlinear state equations 
of motion which are integrated to set the next time 
point. FMAT calculates the partial derivatives of the 
nonlinear state equations with respect to the primary 
states, GMAT calculates partial derivatives with 
respect to the control inputs, GAMAT calculates partial 
derivatives with respect to the process noise variables, 
and BMAT calculates partial derivatives with respect to 
the secondary states (or biases). 

MEAS contains the measurement equations for the appropriate 
estimator. It has entry points YHATV, HMAT, DMAT and 
CMAT. YHATV contains the nonlinear measurement equa- 
tions. HMAT calculates the partial derivatives of the 
nonlinear measurement equations with respect to the pri- 
mary states, DMAT calculates the partial derivatives with 
respect to the control inputs, and CMAT calculates the 
partial derivatives with respect to the secondary state. 

FILLXI is used following either a STATE or MEAS subroutine 
entry point call to create a reduced order (i.e., 
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compressed) vector for use in the subsequent filter/ 
smoother calculations. 

BFLOCS performs the bias-free local smoothing computations. 
Local smoothing iterations are used in conjunction with 
BFFLTR to improve the linearization of the primary states 
in the bias -free extended Kalman filter calculations. 

BFLTR performs the bias filter computations, where the second- 
ary (or bias) state estimates and covariances are 
calculated. 

CONTRL finds the values of the control inputs at any specified 
time by linear interpolation. 

CORECT performs the correction to the bias free state estimate 
and covariance to account for the secondary state esti- 
mates and covariances. This results in the composite 
state estimate and covariance. 

CSLOCS performs the composit state local smoothing computa- 
tions. Backward local smoothing iterations are used in 
conjunction with BFFLTR, BFLTR and CORECT to improve the 
composite state estimate. 

DUMPIT is an output diagnostic routine that dumps the D matrix. 

STORIT is the routine that stores the filter output for subse- 
quent printing and plotting. 

FIS performs the fixed interval smoothing computations. This 
backward smoothing is used on the primary states to pro- 
vide a smoothed state estimate given the entire data 
record. An estimate of the process noise is also calcu- 
lated . 

OUTPUT is the subroutine that controls all program output - 
the printed time history tabulations, printer plots and 
writing data to a tape or permanent disc. Calls subrou- 
tines WRTDAT , PRTDAT and PLTDAT. 

WRTDAT writes time history data to a permanent mass storage 
file. 

PRTDAT prints time history tabulations. 

PLTDAT generates printer plots of the time histories. 
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2.2 KALMAN FILTER EQUATIONS 


The DEKFIS program utilizes a Friedland-Duffy extended ; lman 
filter with locally iterated smoothing to provide filter estimates 
and a fixed interval smoothing algorithm to provide gusts and 
smoothed estimates. The system equations can be written as: 


x x (t) - f(x lf x 2 , u, t) ♦ r(t)w(t) , x 1 (0) • x 1 (1) 

o 

x 2 (t) - 0 , x 2 (0) - b 0 (2) 

y(t) » hCx x> x 2 , u, t) ♦ v(t) (3) 

E{w(t) } « v(t) E{ [w(t) -w(t) ] [w(t) -w(t) ] T } ■ Q(t)fiCt-T) (4) 

E{v(t)} « V(t) E{ [v(t)-v(t)] [v(t)-v(t)] T } » R ( t ) 6 (t-x) (5) 


where x^ is the primary state and x, is the secondary state. The 
nonlinear f and h equations are linearized at each time point (as 
per the extended Kalman filter), defining the following matrixes. 

F ■ (n^ * ) $ - exp {FA t) H * (m * n^ 

8 “ Tlh (n l * n 2^ C " 377 * n 2^ 

r D - f' 1 (j-i) r (n x -< P) 

Using the Friedland-Duffy approach, the filter equations can be 
written as sets of bi^s-free and bias filter equations. 

Bias-Free Filter Equations: 


(1) 

M. » 

X 

t«p » T * r D q r^j., 
/ i . 

(2) 

- 

X 

/ f C x i-i* b i-r u i-i- * • x i-i 

t i-: 
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(3) 


W. - (HMH T ♦ R) t 

(4) K^. - (MH T W* 1 ) . 

■>1 i 

(5) Pi • II - K x H].M. 


(6) 

u i * n • 

hC *i. 0, , u it t i ) 

(7) 

■ *i * 

K Xi V i 

where 

x ■ x^ and 

b • x. 

Filter 

Equations : 


CD 

u i ■ *i-i 

-1 

V i-1 * F i-1 ! * i - 1 ' 1 | B i 

( 2 ) 

s i ■ H i“i 

+ C i 

(3) 

V i - V. - 

K x t Si 

(4) 

w. ■ w. ♦ 
1 1 

S i P b S I 

1 b i-l 1 

(5) 

K h * P\ 

b i u i- 

i S I V 

(6) 

p bi * Cl * 

%i »ll Vl 

(7) 

A /\ 

b i - b i-l 

* K bi ("i • S i b i-l> 


Note that the primary state estimate at the i^h time point is 
denoted by and the secondary state (or bias) estimate is denoted 
by b^. The bias filter also utilizes Kxi» Wi and Vi which are 
computed in the bias-free filter. These two sets are fully coupled 
by the composite state update which corrects the bias-free state 
estimate for the effects of the bias. 


Composite State and Covariance Update 


CD 


x c 

i 


~ 1 

►—4 

< 

i 


x > 

1 



0 I 


. 6 i. 
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( 2 ) 




i 


( 3 ) 


P i * V i P b< 


( A ) 



V i 



For the local smoothing options, it is possible to iterate on either 
the bias-free state estimates and/or the composite state filter 
estimates . 


Bias-Free Local Smoothing: 


d) ViU ' *1-1 * h-i *I-i t 1 - K /lI H I R 'S 


Composite State Local Smoothing: 


c C T I 


*i-m ■ *1-1 * p i-l '1-1 i 1 


fK 1 

) T 

r H T ~ 

X 

_'*V 

[H;C] j 

ct 

i- -J 


R' 1 v 


X 


whert. 


( 3 ) 


$ F* 1 {$- I }B 


0 I 


The fixed interval smoothing equations are used when the filter 
estimates are not adequate. These equations are necessary when 
s ne spates are driven by process noise (i.e., gust states) and it 
is necessary to estimate this process noise. A further benefit of 
fixed interval smoothing is that smoothing estimates are computed 
for all the primary states. These equations are shown on the 
following page. 
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Fixed Interval Smoothing 


(1) x i N * x i ’ P i*. X i 


(2) w. 


i N ' w i * 


T 

l" i i 


(3) A. 


-I * (r • p i H i R i lH i ) M ^ - h.x.,j, 


2.3 STATE AND MEASUREMENT EQUATIONS 

The DEXFIS program has the capability of handling any non- 
linear state and measurement equations that are programmed in 
STATE and MEAS subroutines. The rotorcraft version of DEXFIS ha r 
the equations for three independent estimators: the fuselage/gust 

estimator, the rotor state estimator and the RSRA estimator. These 
equations are discussed below. 

2.3.1 Fuselage/Gust Estimator Equations 

The uncompressed state and measurement vectors for the 
fuselage/gust estimator are shown in Table 2.1. The user has the 
optioi of choosing which states (and hence which equations are to 
be integrated) , which measurements and which process noise sources 
are to be used for a particular problem. 


13 








The X} vector is of size 25x1 and contains the primary states 
for the fuselage/gust estimator. These are: the aircraft 

attitudes - Euler 4,0, <i>; the aircraft inertial velocities in 
the aircraft body axis - v , v , v ; the fuselage angular 

A 7 £ 

rates— p, q, r; the linear accelerometer indicated acceleration - 
a x ,a yi a Zl (this differs from the actual linear acceleration by 
a £irst order lag); the angular acce^ rometer indicated acceleration 
~ ^I*^I**I ( li * ewise lagged) ; gust velocities in an inertial north, 
east, vertical frame - Vg^, v g E * v g v ; and the aircraft position in 
an inertial north, east, vertical frame - Xjj,Y E ,Zy. 

The v vector is of size 24x1 and contains the measurements 

for the fuselage/gust estimator. These are: the aircraft 

attitudes - p ,0 ,<|> ; the aircraft airspeed in the aircraft body 
m m m 

axis - v x m » v *y n »V 2 m (or either alternately or redundantly as 

VmfSm.am measurements); the aircraft angular rates -p,«,,q m »r : 

the linear accelerometer measurements - a XT , a VT , a. T ; the 

*Im yi m “Im 

angular accelerometer measurements - pj m , Al m » fl m I and radar 
position measurements in an inertial north, east, vertical frame 
“ Xn V Y Em’ : V n (or alt ernately V 0R ra . *R m > the radar range, 
elevation and azimuth coordinates.). The associated units are 
to the right of this vector. (The actual measurement data input 
to the program should have these units.) 

The x, vector is of size 50x1 and contains the biases 
(ie b's) and scale factors (ie k’s) for each of the measurements. 

The w vector is a 9x1 process noise vector. This process 
noise drives the linear accelerations - w^w^Wj; the angular 
accelerations - w 4 ,Wg,w 6 ; and the gusts -w 7 ,wg,w g . 

The u vector is a 4x1 control vector that is dummied up 
since there are no control inputs used in the fuselage/gust state 
equations . 

The nonlinear fuselage/gust equations are presented below: 


IS 



Primary State Equations 


1 

I 


♦ • p ♦ (r cos $ ♦ q sin ♦) tan 9 
i • q cos ♦ - r sin ♦ 

1 • (r cos t * q sin #)/cos 9 

v » v r - v q-g sin 9 ♦ g sin 9(o)«- + w $ )l yi 

x y • 

♦ Cq 2 '-r 2 )t xi -qpt yi - rpt^ 

v - v p • v *g sin $ cos 9 - g sin ♦(o)cos9(o)-*-w 2 - w 6 lx 2 

y J X 

♦W 4 ^t - p q 4 X2 * (P 2> r 2 ) iy 2 - rqtj 2 

v -v q • v p ♦ j cos dcos9 - g cos $(o) cos0 (o) + w. - -* 4 t 

2 x y ' 3 

+w S*x3 - pr t x - -q r i yj * (p 2 *q 2 H S3 


t> - w 4 

q * " s 

T - W 6 


a X£ 1 

* a l a xi " 

a i W l 

a YI 

• a 2 a yI - 

a 2 W 2 

*si ' 

" a 3 a *I “ 

a 3 W 3 

(Pi) 

• 

- <Vl - 

a 4 W 4 

(q*j) 

- a 5 qj 

°S W 5 

• * 

(rj) 

• 

■ Vi - 

a 6*6 

V «N 

' a ' V *N 

’ a 7 W 7 

V SE 

‘ a 8 V *E 

' Vs 
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Vg v " °9 V *v * °9 w 


X ■ v cos 9 cos ii» ♦ v u sin 9 sin ♦ cos it - v cos 9 sin i* 
n x y y 

* v 2 sin 9 cos tt cos 9 ♦ v z sin <t sin if; 


■ v x cos e 

sin 

♦ v v sin 9 

sin 

$ sin ♦ v v 

COS o COS it 

♦ v, sin 9 

cos 

$ sin it - v- 

s in 

$ cos # 


■ - v x sin 

6 ♦ 

v v cos 9 sin 

* ♦ 

v 2 COS 9 cos 

* 


Secondary State Equations 



Measurement Equations 

»m ' k , ' * *V n l 
9 . ■ k 8 • 9 *V n 2 

»« * k » *♦ *V n 3 

v x m ■ k v x Cv x * v Xg ) * b Vx * n 4 

'Vra • k Vy ' V y * v yg ) * bv v * "5 

v: m • k v . (v. * Vj ) ♦ b v . ♦ n 6 

p » k • p + b_ * n- 

p r p 


1 ? 


% ’ k q ' <> * b , * " S 
r„ • k r • r ♦ b r ♦ n 9 


**1. ■ k a XI (**i ‘ * s>» i ( 0 )) <• b lxi * n ia 


a 


yi 


- K a - g sin $ (0) cos 6(0)) *b a T ♦ n 

m 1 1 * 1 y I 
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a z Im * k a 2I Ca 2 r - g cos * (0)cos6 (0) ) ♦ba^ ♦ n 12 


> T - k* • p T ♦ b* ♦ n, 7 
*m p I 1 Pi 13 


* k p r ' «i * b qi * " 
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r I * k* * r r ♦ b_ ♦ n. c 

% r I I rj 15 


Xb, » k x ’ XN' > b x ♦ n. A 


V Em ■ k V£ • V E • by E • n r 


’m 


♦ b, ♦ n 


■v 
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V • k 
m K v 


• UVV Xg >q C :4’ r t ‘y 4 ) ‘^V v y * i X 4 ‘P 


* tV : *V P V qi k 4 ):, ' /2 * b v * "19 


am 


-1. 


k 3 • tan Uv *v * r l -p i. s )/ t v x *v ♦q • _ s -r.' v „) 

• .7 * > J 5 - 3 

' ! V n :c 


1 1 


IS 


n 



tinXcvV P V ,, *6 ,/( V v y < *V r V« )) 

*V n 21 

R » ■ k R ' X N 2 * Y E 2 * Z V 2)1/2 * b R * “22 
X ‘ t.n' l (-2 v /(X N J -Y E 2 ) l/2 ) * I> 8r * n 23 

X ‘ V t«' l C-VV * b o R * "24 


where , 


V 

t V v “ 

)cos8cos^ 4 ‘(Vg 
N E 

■V 

cos8sini|/- (Vg. 

+V W )sin 0 
w v 

V 

<V V « 

, )sin9sinqcos^* 
N 

(V *N 

♦■V ) cos8s in 9 

N 



* tv » E *V 

sin 9 

sin o 

sin 

If 

♦ <vv 

cos q COS li; 


* (Vg *v ) 

°V V 

cos 3 

sin o 





v, ■ 

“g 

(V *N * 

V 

sin 9 

cos 0 

cos 


*C'V V WN> 

sin $ sin w 


* X E 

* V w 

) s in 
E 

3 cos 

$ sin 

* - V 

) sin o cos 
E 


* (V Rv 

♦ v w 

) cos 

V 

9 cos 






The parameters used in these equations are defined in the 
fuselage/gust estimator parameter list shown in Table 2.2. The 
terms n., i ■ 1,24 represent random measurement noise. All locations 
are defined in the aircraft body axis system -- positive forward, 
right, down from the aircraft center of gravity. 
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Table 2.2 

Fuselage/Gust Estimator- Parameter List 


PARAMETER 

INDEX 

PARAMETER DESCRIPTION 

1-25 

x^ ic's (use uncompressed vector) 

26-75 

x- Ic's (use uncompressed vector) 

76 

lx i 

77 

iyj a x accelerometer location, ft. (relative to c.g.) 

78 

l *l 

79 

£ *2 ) 

80 

ty 2 > a accelerometer location, ft 

81 

1 

82 

£x 3 ) 

83 

iy- ' a accelerometer location, ft 

84 


85 

* x 4 ) 

86 

ly 4 J V probe location, ft 

87 

iZ4 ) 

88 

ix 5 ) 

89 

tyc ' 3 vane location, ft 

90 

*Z5 1 

91 

£x 6 ) 

92 

iy g oi vane location, ft 

93 

lz 6 ‘ 

94 

a l \ 

95 

®2 1 

96 

a 3 \ accelerometer time constants 

97 

a 4 ( 

98 

a 5 

99 

3 6 ' 

100 

°7 ) 

101 

ap ' gust power spectral time constants 

102 

a 9 1 

103 

v “n ) 

104 

Vw E steady gust components (wind) 

105 

vw v ’ 





2.3.2 Rotor State Estimator Equations 


The uncompressed state and measurement vectors for the 
rotor state estimator are shown in Table 2.3. As in the previous 
estimator, the user has the option to flag which states, measure- 
ments and process noise sources he wants to use. 

The x x vector is a 14x1 vector of primary states used in 
the rotor state estimator. These are: rotor coning angle -So; 

rotor longitudinal flapping -61c; rotor lateral flapping -61s; 
simular fixed system coordinates for the rotor lagging motion - 
Co, clc, els; rotor azimuth angle - and the time derivatives 
of each of these - 6o, 8lc, 81s, ;o, ;lc, els, i R . 

The £ vector is a l^i vector of measurements and associated 

units for the rotor state estimator. These are: the individual 

blade flapping angles 3^ m , i * 1, 6; the individual blade lagging 

angles i * 1, 6; and the cosine of rotor azimuth -cos ip v . 

m 

The x, vector is a 26x1 vector of biases and scale factors 

—t 

for each of the aforementioned measurements. 

The w vector is a 7x1 vector of process noise sources each 
driving one of the seven rotor degrees of freedom modelled. 

The state and measurement equations for the rotor state 
estimator are presented below: 

Primary State Equations 

• I *1 ♦ I H 
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Table 2.3 

Rotor Stat* Estimator-Vector Definition (Uncompressed) 






Secondary State Equations 


*2 " £ 

Measurement Equations 


6 fm ‘ k 81 (6 o ' 8 1 <; ' S l s 5in (V*l )5 * >>6/ 

6 > " k 6i (0 o' 8 l c cos (^R^^i) -8i s sin(g» R ^ i ))+bg i * n i 

0 6 m " k B 6 (B o' 8 l c cos ^r^ 6) ' B l s sin Cv R ^6))*b R6 ^n 6 

«1 B " k Ci (; o ‘ ;i c cos ^R +0 l ) ' i: ls sin(ip R >9 1 ))>b i;i ^n 7 

; i m * k S U o ‘ ; l c co*(* R *ei)-Ci s sin(^ R *0 i ))*b ci ^n i ^ 

? 6m * k C6 C; o * «i cos(’4; R *6 6 )-C 1s sinC* R +0 6 )) ♦ b ♦ 


cosU/ R « k * cos(iO B ) ♦ b ♦ n,, 

m cos^ R C0S ^ R 13 


The parameters used in these equations are defined in the 
rotor state estimator parameter list shown in Table 2.4. 

2.3.3 RSRA Estimator Equations 


The uncompressed state and measurement vectors for the 
RSRA estimator are shown in Table 2.5. As before, the user has 
the option to flag the desired states, measurements and process 
noise sources; but, in this case it is not advisable to do so. 



Table 2.4 

Rotor State Estimator — Parameter List 


PARAMETER 

INDEX 

PARAMETER DESCRIPTION 

1-14 

1c 's (uncompressed vector) 

15-40 

Xg 1c' s (uncompressed vector) 

41-76 

not used 

77 

♦l 

78 

h 

79 

$2 Blade flapping measurement 

80 

, phase angles, rad. 

*4 

81 

>5 

82 

H - 

83 

9 i 

84 

e 2 

85 

9, Blade lagging measurement 

86 

phase angles, rad. 
°4 

87 

*5 

88 

»6 
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I 

1 


RSRA Estimator-Vi 



10x1 


Table 2.S 

or Definitions (Uncompressed) 



1 * 

6x1 



ra«j/s«c 2 

rad/stc 2 

rid/*«C 2 

rjr1/s*C 

rid/tK 

!-ad/s«c 






The integrity of the RSRA estimator requires thet all measure- 
ments be available. 

The vector is a 10x1 vector of primary states for the 
RSRA estimator. These states ire: thr rotor forces at the hub 

in the aircraft body axis system - X R , Y r , Z r ; the rotor moments 
at the hub - L Rt M R , N R ; the inertial accelerations of the 
transmission in the body axis system ~*t x > a t » a t . > and the sum 
of the engine and tail rotor drive shaft torques - C . 

The £ vector is a 10x1 vector of measurements consisting 

of: the transmission load cell reactive forces - A 3 E,C,D,E,F: 

transmission accelerations -at ,a t ,a r _ ; and total drive shaft 

xin y m u z m 

torques - Qt . 

The x ? vector is a ZOxl vector of biases and scale factors 
for each of the 10 measurements. 

The u vector is a 6x1 control vector of body angular accel- 
erations - p,q,r and body angular rates >p,q,r that are treated as 
deterministic inputs. These six quantities should be Obtained 
from a prior fuselage/gust estimator run. 

The w vector is a 10x1 vector of process noise sources that 
drives each of the states. 

The state and measurement equations for the RSRA estimator 
are presented below: 


Primary State Equations 


i 








Measurement Equations (bias $ scale factors not shown) 


X R +h 12 Y R 

. ‘ H 13 Z R 

* h 14 l R * » 

IS M 

R * h l«<JR ‘ 

h ll M 

tat x 

: m t a t y - 

h 13Vt 2 

* 1 14^ t xx^ 

- h 14 (I t -I„ 

14 C 22 tyy 

) r p 


f 6 m t a t y 

- h 14«t - 

h 15 I t rT ^ ' 

h l3 

CI t ^t } 
c xx c zz 

q r 


, f 5 M t a t 2 

h lS f 6 m < 

■ a t * h lA I i- X 

■ lo t 22 

h 16 (It yy' 

l t ) 

t xx 

pq 

f c m t a t 
i 5 1 t v 

4 

* "l 






X R * h 22 

: Y R * h 23“R * h 24 L R + 

h 23 

M R * h 26 Q R 

’ h 21 

Vt 

: M t a t y • 

h 23 M t a t z 

- h-M^t P 
24 L XX 

- h„ 

U 

4 x * t zz ^yy 

) r p 


*«Vv 

' h 24 Q t ‘ 

' h 2S I Vyq * 

h 25 

(I ^xx- It zz ) 

q r 


; f sVt, 

h 2S f 6 m t 

: a t x ' h 26 I t 

• 

r 

zz 

‘ h 26 (I t yy * 

l z > 
c xx 

p q 


ft 26 f S M t a t v * h ; 


f. 


Tnmj R 

ur+Fi)“ R 

TTRI T 

£ 1 a., a. 

l tyy 

> yy q 

♦ ( t XX 

( f 1* 1 2l * * 

= (fl_£T 

TFT 

f, 

a r ♦ 6 t 

a t ♦ n. 



M r 3 

R 


m a t 
t u x 


Ii 


f 5 "t 

m*77) 



1 


where , 


- _h _ Y d - 1 Qn ♦ f 2 mj. *• z zz 

I n>rrr R Tnvrry R 1 V xrnrri 

♦ ^ t yy ^xx^PQ ♦ f 5 m t a «. ♦ a. 

ftfr + jy) r w .-w 1 4 


rmrj v y 


E • fl+f3-£2 Y n - 

"TTT73 R ITT 


1 Q„ - £l + £3-f 2 m^a, - r 

tttt r _ t c y ' ( nrn ) 


* (I ?Y V ' l txx ) pq * t ■ £?- T a t * n 5 
(FI *£3) TTTJT3T y * 


1 


F • X R - m t 3 t ^ • n 6 


a t * a t * ^ 

X m X 

m 


a t " a t * n S 
v V 

• m 


a t. " a t, * n 9 


V * <3 t * n io 

m 


h n ■ h 21 


•1 f 3 

r ttptft 


h i: 


-h,, * -f. * tW*: 
- -i v;rfT"£5) 


23 



h 13 “ h 23 " 'I 


7 TTHTTf 


h. . ■ -h, . • 1 
!•* 24 — 

* t 


h 15 * h 2S 


TTFFnr 


h,. - -h^ - ( f4- f 3 ) 

i5 - 6 y^ rn ^ h T 


The parameters used in these equations are defined in Table 
2.6 (RSRA State Estimator - Parameter List) and in Figure 2.3. 
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Table 2.6 

RSRA State Estimator — Parameter List 


PARAMETER PARAMETER 
INDEX DESCRIPTION 

1-10 x, ic's (uncompressed vector) 

11-30 x r Ic's (uncompressed vector) 

31 y t 

“ f . 

33 f 2 

34 Defined below 

35 f 4 

36 f 5 

37 f 6 

38 Transmission mass, slugs 

39 I. Transmission principle moments of inertia, 

wx slug-ft: 

40 ‘tyy 

41 'm 

• X 

f l * T cos \ * (z R + z t ] sin \ 
x t 

f 2 * T cos \ ‘ ^ Z R + z t^ sin i t 


f 3 * (z R + z t ) cos i t + y- Sin i t 

x t 

f 4 » (z R + z t ) cos i t - y- sin i t 

f s ■ z n sin \ 

U * z r 005 \ 
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Table 2.6 (Concluded) 


and, 

x t* z t "* are the rotor counting geometry dimensions in the 
1 1 1 transmission principal axis system, ft. 

z„ - is tne a « stance from the rotor hub to the trans- 

mission center of gravity *leng the shaf* ft. 

i t - transmission incidence mill, respect to the longi- 

1 tudinal body axis. 









III. PROGRAM DESCRIPTION 


Two forms of input are used by the DEKFIS program. Information 
about the program options, parameter values, initial conditions and 
covariances are read from cards. The time histories of the measure- 
ments and controls are read from a mass storage device such ?.s 
magnetic tape or disk. 

3.1 CARD INPUT 

All data cards are read by subroutine INPUT. The TEKFIS 
input card deck setup is discussed in Table 3.1. This table lists 
the input cards in order according to card type, indicating the 
numbers of each type of card required. Format and variable 
descriptions are discussed for each card type. 

3.2 EXAMPLE INPUT DECK 

An example input deck is shown in Figure 3.1. This deck is 
setup to process CH- 53A flight data through the fuselage/gust state 
estimator. It is to be run in the (standard) combined primary 
and secondary state filter mode without any locally iterated 
smoothing. The output for this deck setup is discussed in Chapter IV. 

3.3 TAPE AND DISK FILE REQUIREMENTS 

The tape and disk file requirements are shown in Table 3.2. 

There is one input file which is required for any filter run. 

This input file contains the measurement (j^) and control (u) 
time histories. There are two output files which are optional . 

These contain state, measurement and process noise estimates 
from a filter or fixed interval smoother run. A define file is 
required for any fixed interval smoother (FIS) run. This file is 
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Tab 1 e 3.1 

DI-KFIS Input Card Setup 
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Table 2.1 (Continued) 
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Tabic 3.1 (Concluded) 


I 

1 

I 


DESCRIPTION 

Refer to card type 12. This set of cards defines the 
Initial bias (secondary states x 2 covariance Matrix - 

P X2 (0) (symmetric matrix). 

Refer to card type 12. This set of cards defines the 
process noise covariance matrix, Q (symmetric matrix). 

Refer to card type 12. This s.‘. of cards defines the 
measurement noise covariance, R (symmetric matrix). 

Number of data channels on the Input tape (default unit 2) 
for a given time point 

M=MY+KU, the total number of measurements + controls 
needed by the program. IDUM is a sorting vector Indicat- 
ing what tape channel is used for the appropriate meas- 
urement or control, i.e., I0UM(3) = 25 Indicates the 25th 
channel on the input tape is used for the third measure- 
ment. 
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figure 3.1 Sample Input Listing 
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Figure 3.1 (Continued) 



DATA CAROS 


1 

1 


ooo 


^ ^ ^ ooo -*-*•*-*— o o o O O m « « O o 

O « N N O C <3 oop OOOOO NNN poooo— ooo^oo 

o o e o «* o o o ooooooooo oooooo o o o o o o — -• o o o ** o o 

oooooooo ooooooooo oooooo oo ooooooooo*»oo 


o 




♦ 


* 


4 


iiguie 3.1 (Concluded) 



1 
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Table 3.2 

Tape and Disk File Requirements 


mm 

mm 

OEFAULT 

LOGICAL 

UNIT 

DESCRIPTION 

COMMENTS 

INPUT 

2 

• Necessary for any filter run 

• Contains measurement ^ and 

control u time histories 

• Usually flight data 

• Written in binary 

• Read in by INREAD 

OUTPUT 

(#1) 

3 

• Optional output 

• Contains time t, measurements 
v , and measurement estimates 

£ 


OUTPUT 

(•2) 

4 

• Optional output 

• Contains up to 4 vector time 
histories 

• Choice of primary state esti- 
mates Xj, accelerations (pre- 
diction cycle) x^, secondary 

state estimates x~ and process 
noise w. 

• Written by WRTDAT 

• WRTDAT is called 
once for each vector 
to be stored and an 
EOF written after 
each 

DEFINE 

FILE 

8 

• Used to store matrices from fil- 
ter runs for use in the fixed 
interval smoother runs 

• Used as incut to FIS run 

• File read/write com- 
patibility handled 
automatically 


I 


3 

















used to store matrixes at every time point in the forward filter 
pass for use in the subsequent backward FIS run. 

Prior to the initial OEKFIS runs, an input file must be 
set up on mass storage. This file must be compatible with the 
FORTRAN code in the INRF.AD subroutine. This input file is to be 
written in binary w’th one record written for each time point. 

Each record contains the measurement and control values at the 
given time point. The FORTRAN code below illustrates this further: 

DIMENSION RECORD (NITEMS) 

DO 10 J « 1, MDP 
WRITE (2) RECORD 
10 CONTINUE 

where . 

MDP - number of data (or time) points 

NITEMS * number of measurement and control channels in 
each record 

RECORD * an array of length NITEMS containing the 

measurement and control values needed by the 
DEKFIS run 

Note: NITEMS can be greater than actually needed fcr the run. 

When running DEKFIS, card type 16 tells INREAD the value 
of NITEMS and which of these NITEMS channels are actually 
to be used. Rearranging is also possible with card type 16. 


44 


/ 



3.4 JOB SUBMISSION ON CDC-7600 MACHINES 


The DEKFIS program was developed on a CDC-7600 and requires 
103777 (octal) small core memory (SCM) and 141S20 (octal) large 
core memory (LCM). 

The JCL (job control language) for a typical DEKFIS run is 
s:*ovn in Figure 3.2. This figure shows ar update, compile, load 
and go sequence. The output files are cataloged after a success- 
ful execution. 
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IV. PROGRAM OUTPUT 


Most of the program's output is written by the printer, 
including optional pointer plots. The program also has the option 
of writing data on mass storage devices (see Figure 3.2) , which 
gives the user the information necessary to make off-line, pen- 
and-ink plots (e.g., Calcomp plots). 

4 . 1 PRINTED OUTPUT 


Examples of program output are shown in Figures 4.1 to 4.7. 
These are explained below: 


Figure 4.1 Contains the program title block and a summary of the 
major inputs specified via card input. These include: 
the mode of operation, model type and size, sample 
step size and number, diagnostic printout flags, and 
variable selection flag vectors. 


Figure 4.2 Contains listings of the primary states, secondary 
states, process noises and measurements used in the 
run. A summary table of the output increments is also 
printed out. The remaining information is used primarily 
for diagnostic purposes. 

Figure 4.3 Shows the parameter vector as input on card type 11, 
the initial covariance of the primary states (PS) as 
input on card type 12, and the initial covariance of 
the secondary states (PB) as input on card type B. 


Figure 4.4 Shows the covariances of the process noise and measure- 
ment noise as input by card types 14 and 15. It also 
contains printout generated by INREAP as it reads the 
input data file into local storage. 
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Figure 4.3 (Continued) 
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Figure 4.4 Process Noise and Measurement Covariance 



Figure 4.5 Contains printout from INREAD. The initial time 

points of the measurements and controls used in the 
run are printed out in order to allow verification 
that the correct data is bain.; used. 


Figure 4.6 Contains the intermediate printout of the performing 

filter. This is level 2 printout containing the state 
and bias estimates and the associated covariances at 
the indicated iteration. The level 3 printout (not 
shown) is extremely detailed giving all the inter- 
mediate steps in each iteration. It is not recommended 
that level 3 printout be used except for diagnostic 
purposes . 

Figure 4.7 This figure shows a small portion of the printer plots 
generated by this run. The variables are in groups 
of four and are plotted at the rate of one line of 
printout for each time point as dictated by the print 
increment . 


4.2 MASS STORAGE OUTPUT 


The two types of mass storage output are shown in Figure 3.2. 

It is important to be able to read these files in order to set up 
runs for plotting, etc. 

Output File »1 - This file contains time, the measurements, 
and the measurement estimates. The program control logic for 
writing this file is conveyed to the DEKFIS run on card types 9 
and 10. The file is written by subroutine OUTPUT, and the following 
FORTRAN code is supplied for reading the file: 
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DIMENSION TITLE (18) , YHAT(MY), YM(MY) , Ui‘(MU) 
READ (3) TITLE, MY, MU 
DO 10 J-l, NDPF 
READ (3) T, YHAT , YM, UU 
1G CONTINUE 


TITLE 

- is 

the 

title information 

YHAT 

- is 

the 

£ vector of length MY 

YM 

- is 

the 

^ vector of length MY 

UU 

- is 

the 

u vector of length MU 

MY 

-is 

the 

number of measurements 

MU 

- is 

the 

number of controls 

NDPF 

- is 

the 

number of data points on file 

T 

- is 

time 


Output File *2 - This file can contain up to four vector time 
history files, written sequentially and separated by end-of-file 
marks. The user has the option of specifying which (if any) of the 
following four vector time histories are to be stored. These are 

A 

(and in the order stored): the primary state estimates , the 

primary state accelerations (after the predict cycle of the Kalman 

A 

filter iteration) x^, the secondary (or bias) state estimates 
and the process noise estimate w. The vector time histories are 
each written by WRTDAT as specified b> the options of cara 
typo 10. An end-of-file mark is written by WRTDAT, so an EOF 
separates each section. The following FORTRAN is supplied for 
reading any one section: 

DIMENSION TITLE(18) , VECTOR (M) 

READ (4) TITLE, M 
DO 10 J * 1, NDPF 
READ (4) K, VECTOR 
10 CONTINUE 
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Figure 4.5 INRFAD Output 
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l ; igure 4.6 (Concluded) 
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where. 


TITLE - is the title information 

A A 


VECTOR - is 
M - is 
NDPF - is 
K - is 


a vector of either 
the length of the vector 
the number of data points 
a dummy integer variable 


, £2 or - 

on each file section 
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V. EFFECTIVE USE OF THE DEKFIS PROGRAM 


The pexformance of the DEKFIS program is dependent on the 
proper choice of the initial conditions for the primary and 
secondary states (x and b_) , the initial covariances associated 

— O "V 

with the primary and secondary state initial condition estimates 
(P and ) , and the process and measurement noise covariance 
matrixes (£ and R) . The following systematic procedure has been 
developed for determining these program inputs: 

1 . Compute measurement noise covariance matrix-R. Assume 
R is diagonal and calculate for each measurement. This 
can be done approximately by inspecting the measurement 
time history and picking out the peak-to-peak value of 
the zero mean gaussian measurement noise on that channel. 
Assume this peak-to-peak value = 60 

2. Estimate values of biases and scale factors for each of 
the measurements. Sometimes these are readily apparent. 
If they are not known, set the biases equal to 0.0 and 
the scale factors equal to 1.0. 

3. Estimate the initial covariance (Pfc ) for each of the 
biases and scale factors to be identified (i.e., treated 
as secondary states). Assume Pk is diagonal and 
calculate each variance term (a*j , keeping in mind that 
6 is the total peak-to-peak excursion allowable. 

4. Estimate the initial conditions on the states. Using a 
best estimate for the measurement at t»0 (usi lly the 
initial data point is as good as any) and the oias and 
scale factor estimates determined in step 2, back 
calculate through the measurement equations to estimate 
the states at t»0. 

5. Estimate the initial covariance (P Xo ) f° r the initial 
conditions on the states. This is possible by taking 
the measurement noise covariance from step 1 and covari- 
ances for the biases and scale factors from step 3 and 
using the measurement equations for back calculation. 

6. Estimate the process noise covariance matrix - 0. Assume 
2 is diagonal. The estimate for £ is the most difficult 
to determine, and several runs may be required to give 
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satisfactory results. A low value for a given diagonal 
element of £ will result in a false convergence of the 
states, and~the measurement matches will show discrep- 
ancies. A large value for a given element of £ will 
lead to unconverged states, and the update correction 
in the predict/update cycle of the Kalman filter will 
be large. In other words, the measurements will drive 
the filter leading to deceivingly good measurement 
matches. Ir is necessary to find a value for 0 which 
will give good measurement matches, but also wrll be 
converged. 

This procedure was used in the validation of the rotor state 
estimator. One example case was run assuming only blade flapping 
degrees of freedom. Lag degrees of freedom were omitted by simply 
setting input flags. Thus, the state variables were: 

- T " t 8 o ,8 lc ,3 ls , '' , 'r ,8 o ,8 lc’ 8 ls’’ 4 'r^ 


The example was for a three-blade rotor, and the measurement 
variables consisted of the flap angle of each blade and the 
cosine of the rotor azimuth angle; i.e., 


y i" 8 i* B o‘ 6 lc cos( ' 4 'r +l,, i ) ' 8 ls sin ^r +<,, i ) * 

i * 1,2,3 


y 4 - cos O r ) 


where the blade spacing angles, are known. The measurements 

were corrupted by relatively small amounts of white noise. 

The inputs to the filter were the state initial conditions, 
x Q , the initial state covariance matrix, p Q , the measurement 
noise covariance, R, and the process noise covariance, Q. 

When processing actual flight data, x 0 , P Q and R are reason- 
ably well known. Because the performance of the filter largely 
depends on the ratio of the norms of P , R, and Q, it is the 



choice of the less -well* known Q which determines the performance. 
The rotor estimator is simply: 


1 

X* 
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where 


E [w(t.)] - 0, E[w(t i ) w( tj ) T ] - Q6 ij 
and (in this case) 

-1 ^ 3 o’ 8 lc’ 3 ls’ ‘ 


It can be seen that setting Q to zero tells the filter that 
the states are constant and that letting Q approach infinity 
(relative to R) tells the filter to weight the measurements very 
heavily. When estimating trim tip-path-plane angles and instru- 
mentation biases, Q may be set to zero. When the tip-path 
plane is not constant due to a control input, an intermediate value 
of Q is desirable. Too large a value for Q results in the state 
estimates being heavily driven by the measurements. 

The effect of Q on the rotor state estimator was investi- 
gated for this example problem. Figure 5.1 shows a portion of 
the measured and estimated time histories of the flapping angle 
of the first blade from a run for which Q was set to zero. A 
control input beginning at t»2.S seconds causes a large change 
in the blade’s flapping response, but this change is not tracked 

by the filter. Figure 5.2 is the same plot from a run for which 

0 

Q was set to 10 . Here, as expected, the measurement was tracked 
exactly (the symbol ”B" overwrote the symbol "A" everywhere) . 
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1 
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I'igure 5.2 Rotor State Estimator Example with Q « 10 
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Figure 5.3 is the same plot from a run for which Q was set to 
1.0, and it can be seen that the measurement is partially tracked. 
A somewhat higher value of Q would be required for satisfactory 
results . 
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